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Abstract 

We describe the simulation data produced by a pilot programme to compute mock 
weak gravitational lensing maps for a range of currently popular cosmological mod- 
els by ray tracing through high-resolution N-body simulations. The programme 
required only a modest investment in computer time to produce maps accurate 
to arcminute scales covering hundreds of square degrees of sky for 4 cosmological 
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[ 1 Introduction 

Weak gravitational lensing by large-scale structure has become an indispens- 
^ I able tool for modern cosmology, used already to set constraints on the mass 

(73 ■ density {flmat, in units of the critical density) and the fluctuation amplitude 

(cTg) (see e.g. van Waerbeke & Mellier, 2003; Hoekstra, Yee, & Gladders, 2002, 
for the current status) and touted for its potential to constrain cluster scaling 
relations (Huterer & White, 2002) and dark energy (Benabed & Bernardeau, 
2001; Huterer, 2002; Hu, 2002; Heavens, 2003; Abazajian & Dodelson, 2003; 
Refregier, 2003; Jain & Taylor, 2003; Bernstein & Jain, 2004; Takada & Jain, 
2004; Takada & White, 2004). Like the anisotropics in the cosmic microwave 
background (CMB), the theory of weak gravitational lensing is well understood 
(Mellier, 1999; Bartelmann & Schneider, 2001). Unlike the CMB however the 
calculation involves modeling the non-linear evolution of the mass in the uni- 
verse. This makes the predictions of weak lensing very rich, but also means 
that an accurate treatment requires numerical simulations. 
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In this paper we give some details of a pilot program designed to provide weak 
lensing ray-tracing simulations of a small grid of cosmological models in the 
currently favored region of parameter space. These models can be used by 
observers wishing to test their algorithms or fit to their data and by theorists 
wishing to test or calibrate fast and flexible, but approximate, methods of 
calculation. The initial grid is small, but the process of its creation is almost 
entirely automatic allowing it to be expanded simply as the need arises. 

This paper focuses on the creation of the model grid. We describe the choice 
of cosmological parameters in §2.1, the N-body simulations in §2.2 and the 
ray tracing and numerical issues in §2.3 and §2.4. We give some representative 
results in §3 before concluding in §4. 



2 Cosmological models 

2. 1 Choosing parameters 

Wc choose a small number of cosmological models which provide good fits 
to the current CMB and large-scale structure data. For simplicity all of the 
models in the pilot program are variants of the cold dark matter model with 
a dark energy component, and their parameters are shown in Table 1. The 
parameter variations around a base model are designed to keep the CMB 
fluctuations almost invariant, in anticipation of increasingly precise data from 
WMAP and Planck, and the entire process is automated using Perl scripts 
and C code. 

We begin with a relatively small number of models to demonstrate the cost 
and feasibility of making such grids. As we gain experience in using these 
data products and understand the drivers better we can extend the model 
grid and/or increase the fidelity of the simulations. Are our priorities to have 
larger maps? higher resolution? more redshift range? more 'sky' ? etc. 

We first pick the physical matter density Wmat = ^mat^^, the (comoving) 
distance to last scattering, rfis, and a (constant) equation of state of the dark 
energy, w. We approximate the redshift of last scattering as z = 1080, ignoring 
the slight matter and baryon density dependence. This then allows us to solve 
for the Hubble constant, /i, and thus ilmat and Qde assuming spatial flatness. 
Because they are reasonably well known, compared to some other parameters, 
we fix cUmat = 0.145 and d\s = 13.7 Gpc for all of the models in our grid. These 
numbers are close to the best fit for a recent analysis of WMAP and SDSS 
data (Tegmark, 2003). 
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Model Ornat w h n as T 



1&2 



0.296 -1.0 0.70 1.00 0.93 0.15 977 



3&4 



0.357 -0.8 0.64 1.00 0.88 0.15 975 



5&6 



0.296 -1.0 0.70 0.95 0.85 0.10 979 



7&8 



0.357 -0.8 0.64 0.95 0.81 0.10 976 



Table 1 



Parameters for the models run. For each cosmological model two independent sets of 
initial conditions are generated. All models are spatially flat, so fl^e = 1 ~ ^mat- For 
all models the matter density is tOmut = 0.145 and the baryon density ujh = Q.bh?' was 
fixed at 0.023. Our results are very insensitive to this latter choice. All models have 
power-law spectra (no running) with index n and the dark energy has a constant 
equation of state w. The normalization, cts, is from a fit to the WMAP TT power 
spectrum data and the is for this fit with 893 degrees of freedom. 



The next step is to specify the other model parameters, for example the optical 
depth to Thomson scattering, r, the spectral index, n and the baryon density 
uj}j. (The optical depth enters primarily through its effect on the normalization 
of the power spectrum when we fit to WMAP.) We currently hold the number 
of light neutrinos fixed, and include no massive neutrinos. We deal with pure 
power-law spectra with no running spectral index. Our primary variations are 
in the optical depth, the spectral index and the equation of state of the dark 
energy (see Table 1). The first two affect the amphtude and shape of the 
density or potential fluctuations at late times while w is a parameter of great 
interest to the cosmology community. We use a fixed equation of state in this 
initial survey, though it is easy to include any known functional form in the 
future. 

For a specified model we compute the CMB anisotropy spectrum using v4.5 
of CMBFAST (Seljak & Zaldarriaga, 1996). The likelihood software provided 
by the WMAP team (Verde et al, 2003; Hinshaw et al., 2003; Kogut et al., 
2003) and the known anisotropy spectrum are then used to find the best fitting 
normalization of the power spectrum, as, which is then converted into internal 
code units for the N-body simulations. In this way a file of model parameters, 
shown in Table 1, is built up. The lensing convergence angular power spectra 
for our 4 models, assuming the Born and Limber approximations, are shown 
in Fig. 1. Note that the spectra for each pair of models have similar shapes, 
differing primarily in amplitude. In each case the model with the most distant 
sources has the highest amplitude, even though it has a slightly lower matter 
power spectrum amplitude, erg. There will also be differences in the higher 
order moments and the source redshift dependence of the power spectra. 
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Fig. 1. The angular power spectra predicted for our 4 cosmological models for sources 
at Zg = 1. We assumed the Born and Limber approximations, and used an ansatz 
for the non-linear mass power spectrum to make these predictions. The dotted lines 
rising steeply to top right are the shot-noise power spectra for 7rins = 0-4 and 
flgai = 25, 50 and lOOarcmin"^. 

2.2 N-body simulations 

Another Perl program reads this cosmological parameter file and generates 
a directory structure containing the requisite auxiliary files for the N-body 
simulation, plus scripts for generating the initial conditions and running the 
N-body code on the IBM-SP Seaborg at NERSC. 

All simulations require choices to be made about box size, resolution and sam- 
pling. To minimize numerical artifacts we want to run the largest box with 
the highest dynamic range consistent with the computational resources. Since 
this is a pilot program we have chosen simulations which use 256^ particles 
in a 200 h~^Mpc box. As we shall discuss in §2.4, this is expected to be suffi- 
cient for our purposes now, but will need to be improved upon in the future. 
As further resources become available we can increase the volume simulated 
and/or decrease the force softening. 

The initial conditions are generated by displacing equal mass particles from 
a "fuzzy" Cartesian grid a.t z = 50 using the Zel'dovich approximation. The 
initial fluctuations are Gaussian, however we choose from a larger number of 
initial conditions those whose low few /c-modes are close to the mean. This 
avoids running any simulations which are significant outliers and may skew 
statistics for a small number of runs. For each cosmological model we choose 
two independent initial conditions. Ideally we would run more realizations of 
each model to increase the volume simulated, but with two runs we will be 
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able to simulate several hundred square degrees which is enough for now. 

The TreePM code, described in White (2002), is used to evolve the particles to 
z — with the phase space information output at selected times after z — 3. 
The outputs are spaced equally in conformal time, with a spacing equal to 
the time taken for light to travel 50 h~^Mpc (comoving). This spacing is small 
enough to obtain good lensing maps - a halo with "typical" velocity 300km/s 
will move only 50 /i~^kpc between outputs, and its departure from straight line 
motion will be small. Each simulation gives between 80 and 90 phase space 
dumps, depending on the cosmology, each of ~ 500 MB. The gravitational 
force softening is of a spline form, with a "Plummer equivalent" softening 
length of about 30 /i~^kpc, again small enough for our purposes. 

Since the box size is only 200/i~^Mpc, leading to a relatively coarse samphng 
in k, we generate the initial conditions from a fit to the transfer function which 
has the baryon oscillations "smoothed out" . For now we use the fit published 
by Eiscnstcin & Hu (1999). This fit was made to an older version of CMB- 
FAST, and since that time improvements in the code and the physics (Seljak 
et al., 2003) have changed the predictions for T{k) slightly. The differences are 
small, and for the purposes of testing and calibrating algorithms it is enough 
to know the input T{k) precisely, but in the future we will improve this aspect. 



2.3 Making lens maps 

We make two sets of maps. The first, which we shall call the Born series, 
assumes that the convergence field, k., is simply the integral along the line of 
sight of the density field weighted by a simple kernel 



where S is the overdensity, a is the scale-factor, x is the comoving distance 
and g{x) is the lensing weight 



for sources with distribution p{Xs) normalized to J dp = 1. This set of maps 
has the advantage of not requiring a Fourier transform in its construction, 
providing higher resolution for a given pixelization. Its principal disadvantage 
is the approximate nature of the calculation. 

We make several sets of 16 quasi-independent 3° x 3° maps with 1024^ pixels. 
The first set of maps has sources fixed at = 1, while the second set uses a 




(1) 
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source distribution of the form (Brainerd, Blandford & Smail, 1996) 



dp 
dzs 



oc z^exp -{zs/zq) 



(3) 



for zq — 2/3 and 1. For this distribution {z) — r(|)2;o — 1-5 -Zq, but because 
there is more of the lensing weight at higher z the lensing signal is smaller for 
this distribution than a (^-function distribution with Zg — 1.5 zq. 

The second set of maps comes from a full multi-plane ray tracing algorithm as 
described in Vale & White (2003). We make maps of the shear components, 
7j, and the convergence, at a range of source redshifts from 2; ~ to 3 in 
steps of 50 /i~^Mpc. Again we make 16 maps per simulation, choosing different 
positions and orientations for the 'observer'. In each case a 2048^ grid of rays 
subtending a field of view of 3° is traced through the simulation, with the 
required Fourier transforms being done on a 2048^ grid. All assignments to and 
from the Fourier grid are done with cloud-in-cell (CIC) assignment (Hockney 
& Eastwood, 1980). 

The two shear components and the convergence are output at each source 
plane. From this a map of the (measurable) reduced shear, 7^/(1 — k), for any 
source distribution can be computed. For a distribution dp/dzs the weight 
given to source plane j is 



where Ax = 50 /i~^Mpc is the spacing between outputs. The final shear at each 
point is then 7, = Y^j '^jli''^ with 7^^ the ith component of the shear measured 
on plane j. We show an example of the weighting factors for = 1 in Fig. 2, 
along with the (trivial) weight for a 5-function source distribution at 2;^ = 1.5. 
The integral over Xs is extremely well approximated by a sum over the output 
planes. These maps can thus be combined to produce lensing maps for a wide 
range of source distributions, allowing source tomography to be used. It is 
straightforward to introduce mock 'galaxies', with an appropriate distribution 
of intrinsic ellipticities, at this stage and produce mock catalogues. On each Xs 
slice, the angular positions of the galaxies can be made to trace the dark matter 
with some fidelity, giving rise to source clustering effects. Reasonable choices 
produce angular correlations and redshift distributions for the 'galaxies' which 
match current data well. However the parameter space of such catalogues is 
large, so we have chosen instead to make available simple co-added shear maps 
which can be post-processed with a wide variety of simulated galaxy properties 
to make shear maps. To save space we provide the maps downsampled from 



wj = — H{zj)Ax 




(4) 



2048^ to 10242. 
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Fig. 2. An example of the weights to be apphed to the source planes of model 1 
in order to construct maps appropriate to Eq. (3) with Zq = I (solid) and a trivial 
Zg = 1.5 (dashed). The planes are spaced by Ax = 50/i~^Mpc. The weight for 
Zs ^ 1.5 is unity in a single bin. The area under the solid histogram is also unity, 
but 4% of this area lies beyond z = 3 where our simulations stop. 

2.4 Numerics 



To understand the fidelity of the map making procedure we need to account 
for the finite volume, force and mass resolution of the N-body simulation and 
the effect of the finite samphng during ray tracing. The comoving distance to 

our furthest sources is 3 Gpc, at which distance our box subtends 3.8°. The 
volume of space probed by a lensing field is ~ lO^Mpc^, comparable to the 
volume simulated in each run. 

In Vale & White (2003), we developed an analytic model to estimate the effects 
of finite force and mass resolution and the Fourier grid used in the ray tracing. 
This model was based on modifying the 3D dark matter power spectrum that 
appears in the Limber integral for Cg in a way designed to fit the effects seen 
in a series of N-body simulations. This allows us to compare a given numerical 
configuration to a "perfect" model. Fig. 3 illustrates these effects individually 
for model 1 of Table 1 with sources at Zs = 1. For simplicity we use the 
power spectrum fit of Smith ct al. (2003) in computing these ratios. Note that 
the suppression of power by the finite force resolution and FT grid is partly 
countered by the artificial increase in power due to shot-noise. Based on Fig. 3 
we expect our maps to be accurate to several percent, in the 2-point function, 
for i < 3000 and even higher in some situations where the cancellation is 
accurate. 

While the 16 maps per model are not fully independent, they do sample dif- 
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Fig. 3. Resolution effects in the power spectra of tfie maps, compared to the sample 
variance in a 3° x 3° map. The y-axis shows the ratio of lensing power spectra 
to a fiducial model prediction. All calculations arc done with the 3D mass power 
spectrum fit of Smith et al. (2003), the Limber approximation and the model of 
Vale &; White (2003). The dotted line shows the effect of the finite force resolution 
in the N-body simulation, the short-dashed line adds the finite mass resolution. The 
long-dashed line shows the effect of the FT grid in the lensing algorithm while the 
solid line combines all of these effects. The error bars indicate the sample variance 
expected in a 3° x 3° field. 

ferent regions of the box in different stages of evolution and with different 
projection effects. The fraction of the volume traced by the rays in the field 
of view, weighted by the contribution to the variance of the convergence from 
each box, is between 10-20% (for 100 < I < 3000). This forms a rough upper 
limit to the degree of correlation between the maps. Including the two inde- 
pendent boxes per model we have simulated several hundred square degrees 
of sky per cosmology, which is close to the amount of observational data cur- 
rently existing. Extending this to more sky is a straightforward exercise that 
simply requires more computing resources. 



3 Results 



These simulations and lensing maps can be used for a variety of purposes. 
Here we simply show some preliminary analysis to orient the reader. The real 
uses of these maps will be in testing analysis algorithms and new ideas. 

Before turning to the lensing maps let us consider the underlying mass distri- 
bution. Fig. 4 shows the 3D matter power spectrum from our 8 boxes dX z = 1 
and z — Q compared to the semi-analytic model of Smith et al. (2003). The 
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Fig. 4. A comparison of the 3D mass power spectrum with the fitting function of 
Smith et al. (2003) at selected redshifts. The y-axis shows the ratio of the numerical 
to the semi-analytic result, with the two independent runs per cosmology and two 
different redshifts shown as different symbol types in each panel. Open symbols are 
for z = and starred symbols are for 2; ~ 1. The top panel is models 1 and 2, then 
3 and 4 and so on. The vertical dotted line marks the fundamental mode of the 
simulations. 

agreement is good for both the ACDM models and the models with equation 
of state w = —0.8, as expected. On the scales simulated the dark energy only 
enters through the Hubble parameter in the drag term, slowing the growth of 
large-scale fluctuations when the dark energy comes to dominate the expan- 
sion. Since the fit of Smith et al. (2003) has been tested for open and ACDM 
models it is no surprise that it works equally well for w — —0.8. The dis- 
agreement at the highest k might be due to numerical artifacts in the N-body 
simulation or in the computation of the power spectrum, however these points 
are significantly above both the shot-noise and force smoothing limits. We 
find similar behavior in a number of higher resolution simulations we have an- 
alyzed, and a similar level of disagreement can be found in some of the figures 
of Smith et al. (2003) for CDM models and may reflect inherent inaccuracy in 
the fitting function. A qualitatively similar level of agreement was also found 
between particle-mesh simulations and the semi-analytic model by Ishak et al. 
(2003), though they were unable to probe very high /c-modes due to limited 
force resolution. 

As a further check on the N-body code we investigated the scaling of the 
power spectrum of an n = —1 self-similar model with the same numerical 
parameters as the models described. This suggests that the power spectrum 
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Fig. 5. The angular power spectra predicted for model 1 with = 1. The symbols 
with error bars are the mean and variance of the spectra from 16 maps, each 3° x 3°, 
from runs 1 and 2. The points have been offset slightly (horizontally) for clarity. 
The open and filled squares are the results from our "Born" series while the open 
and filled circles are for the full ray tracing. 




100 1000 10000 

I 

Fig. 6. The angular power spectra predicted for model 3 with = 1. The symbols 
with error bars arc the mean and variance of the spectra from 16 maps, each 3° x 3°, 
from runs 3 and 4. The points have been offset slightly (horizontally) for clarity. 
The open and filled squares are the results from our "Born" series while the open 
and filled circles are for the full ray tracing. 

should be accurate to a few percent on the range of scales of interest. 

The K angular power spectrum produced from 16 maps each for models 1 & 2 
is shown in Fig. 5. For this comparison we show maps made assuming the Born 
approximation and with full ray-tracing, both for sources fixed aX Zs = ^- The 
agreement is generally good, with the variance at low-^ expected from sampling 
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Fig. 7. The skewness of the convergence, k, for models 1-8 for sources at = 1. 
The symbols with error bars are the mean and variance of ^3 from 16 'Born' maps, 
each 3° X 3°. The results from the ray tracing are very similar and are omitted. The 
points have been offset slightly (horizontally) for clarity. 

16 3° X 3° fields. The open circles, corresponding to model 2, are slightly higher 
than the semi-analytic model at intermediate i, as expected from the excess 
power seen in the top panel of Fig. 4 (the triangles) at intermediate k. At the 
highest i the 'Born' maps show a slight excess power which could be due to 
shot-noise in the simulations or to a breakdown in some of the approximations 
we have made. Alternatively the semi-analytic model and the ray tracing could 
be underestimating the power on these scales, as suggested by Figs. 3 and 4. 
In measurements these scales would be significantly affected by finite galaxy 
ellipticitics and densities, so this disagreement is likely unimportant. Fig. 6 
shows the same comparison for models 3 and 4. The power spectra for the 
remaining four models are very similar and we do not show them explicitly 
here. 



The maps are non-Gaussian, as expected, and we provide the skewness 



Qli — 



(5) 



for each of the models in Fig. 7. The k field is smoothed by a boxcar of side 
length 6 before the moments arc computed. The plot shows the average and 
the variance across the 16 maps as before. As noted previously (White & Hu, 
2000; Vale & White, 2003) the skewness suffers from a large sample variance. 
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Fig. 8. The shear 3-point function, C+++) foi' isoceles triangles with Oi = 02 = 2' 
as a function of the cosine of the included angle. Different symbol types represent 
the mean and standard deviation of 32 maps, each 3° x 3°, for the 4 cosmologies 
simulated: models 1 &: 2 (open squares), 3 & 4 (filled triangles), 5 &: 6 (open circles), 
7 & 8 (filled squares) . Points have been offset slightly horizontally for clarity. In each 
case we assume z = \. 




Fig. 9. The shear 3-point function, C+++) for isoceles triangles with 6i= 62 = 2' as 
for Fig. 8, except scaled to unity at maxmimum. Note that the shape of the curves 
is very similar among all of the models. Symbol types are as in Fig. 8. 

We also show the shear 3-point function for isoceles triangles of side length 
Oi = 02 = 2' in Fig. 8. Here we plot 10^C+++ ^ ^ function of the cosine of the 
included angle for the 32 fields, each 3° x 3°, simulated for each cosmology 
assuming the sources are at Zg = 1- The points are the mean of the 32 fields, 
and the error bars represent the variance. The highest 3-point function is 
for models 1 & 2, which have the highest normalization (as) of all of the 
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models. Models 3 & 4 have the next largest amplitude, and the second highest 
normalization. For models 5 & 6 the effect of a higher normalization isn't 
enough to overcome the less negative w of models 7 & 8, which have a slightly 
higher amphtude. Notice the large scatter from field to field, but the familiar 
shape (Takada & Jain, 2002). Apart from an amplitude change the shapes of 
the curves are extremely similar, as shown in Fig. 9. We find a similar behavior 
for other configurations we have checked, which has interesting implications 
for the signal encoded in the non-Gaussianity of the maps and how best to 
extract it. We also note that the peak of C+++ appears to be shifted shghtly 
from the equilateral configuration, due to the asphericity and substructure of 
the dark matter halos (Ho & White, 2003; Dolney, Jain & Takada, 2004). 



4 Conclusions 

Gravitational lensing has made rapid advances in recent years, and observa- 
tions are now pushing the limits of existing theoretical models. If lensing is to 
fulfill its promise as a precision tool, theorists need to improve the predictions 
of cosmological models for the n-point correlation functions. An important 
step in this program is the creation of grids of N-body models of sufficient 
resolution and sampling to enable accurate simulation of weak lensing. This 
paper describes a first step in this direction. 

Based on analytic arguments developed in Vale & White (2003) we believe that 
simulations with 256^ particles in boxes of side 200 h~^Mpc are sufficient to 
produce lensing maps accurate to ^ ~ 2000—3000 or scales of a few arcminutes. 
To achieve good convergence in the multi-plane ray-tracing algorithm that we 
use requires time dumps spaced equally in conformal time with a spacing close 
to SO/i^^Mpc. 

The total computational cost of this project was very modest. The major 
human cost is the time spent managing the disk space, archiving and retrieving 
the data to ensure that the total usage stays below quota. The simulations 
each took 2000-3000 time steps, requiring slightly under one day each on 
32 processors of the IBM-SP Seaborg at NERSC, i.e. between 640 and 720 
CPU hours. The phase space data required ~ 40 GB per model while the 
ray tracing simulations required around 100 GB of intermediate storage and 
100 CPU hours per model. The disk usage is driven by the need for all the 
phase space data and the 80 — 90 high resolution source planes for each of 
16 maps. The final stacked and downsampled maps are < 1 GB per model 
and source distribution. The entire process was monitored and controlled by 
scripts, making it easy to extend the grid as more resources become available. 

We have made the raw maps, along with some auxiliary data products, freely 
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available to the community at littp://mwliite.berkeley.edu/ in the hope that 
they will be useful in taking the next step. 

M.W. thanks Tzu-Ching Chang, Dragan Huterer, Bhuvnesh Jain, Jason Rhodes 
and Masahiro Takada for helpful comments on an earher draft. The simulations 
used here were performed on the IBM-SP at the National Energy Research 
Scientific Computing Center. This research was supported by the NSF and 
NASA. 
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